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Preface 


This  report  documents  the  underlying  equations  and  methods  of  the 
DAY-24  model  that  is  being  developed  by  the  U.S.  Army  Research 
Laboratory,  Survivability  Lethahty  Analysis  Directorate  (SLAD)  to 
simulate  infrared  signatures  of  natural  background  scenes  as  affected 
by  atmospheric  and  meteorological  processes.  The  work  is  funded  as 
part  of  the  SLAD  mission  in  support  of  the  Tools,  Techniques,  and 
Methods  thrust  in  infrared  synthetic  scene  generation  and  was 
performed  entirely  in-house.  The  model  was  developed  to«meet  the 
requirements  of  SLAD  System  Leaders  in  tasks  to  evaluate  various 
countermeasure  technologies  being  developed  for  hit  avoidance  and 
target  acquisition. 


iii 


Preceding  Pagers  Blank 


Acknowledgements 


We  acknowledge  the  contributions  of  the  Survivability/ Lethality 
Analysis  Directorate  Field  Team  for  conducting  the  field  experiments 
and  collecting  the  data  for  model  verification.  Specific  individuals  to 
be  acknowledged  are  Lon  Anderson,  Team  Leader,  Steve  Lacy,  Met 
Data,  and  Thelma  Chenault,  Data  Analysis.  Others  include  Robin 
Lacy,  Rebecca  Phillips,  Mark  Nelson,  and  Joe  Trammel,  all  students  at 
New  Mexico  State  University.  Scarlett  Ayres  served  as  the  liaison  for 
contributing  technical  details  from  the  U.S.  Army  Corps  of  Engineers 
Smart  Weapons  OperabiUty  Environment  Program  and  provided  the 
baseline  thermal  scenes. 


V 


Contents 


Preface 

iii 

Acknowledgements 

V 

Executive  Summary 

ix 

1.  Introduction 

10 

2.  Earth-Air  Surface  Formulation 

5 

3.  Equilibrium  Surface  Flxrxes 

9 

4.  Equilibrium  Surface  Temperature 

15 

5.  Time  Dependence 

17 

6.  Application  To  IR  Scene  Generation 

19 

7.  Summary 

23 

8.  References 

25 

Appendix 

A.  Parameter  Definitions . 

. 27 

Distribution 

31 

vii 

Figures 


1.  Sketch  of  boundary  layer  model . 2 

2.  Expanded  sketch  of  earth-air  interface  region . 5 

3.  DAY-24  meteorological  input  and  modeled  fluxes . 10 

4.  DAY-24  effect  of  soil  thermal  properties . 12 

5.  IR  scene  baseline  temperature  map  (upper)  and  soil  map 

(lower) . 19 

6.  IR  scene  emissivity  map  (upper)  and  reflectivity  map  (lower) . 20 

7.  IR  scene  diurnal  cycle . 21 


Executive  Summary 


This  report  documents  the  technical  development  and  validation  of 
the  U.S.  Army  Research  Laboratory  (ARL),  Survivability/ Lethality 
Analysis  Directorate  (SLAD)  thermal  infrared  (IR)  scene  model, 
DAY-24.  The  model  was  developed  in  response  to  a  need  in  SLAD  to 
predict  electromagnetic  IR  signatures  of  natural  backgrounds  of 
various  types  under  various  meteorological  and  atmospheric 
conditions.  The  output  of  DAY-24  is  (1)  a  detailed  "stand  alone" 
surface  temperature  and  "energy  balance"  single  point  analysis,  and 
(2)  an  extended  area  IR  scene,  or  "signature,"  which  can  then  serve  as 
a  baseline  and  input  for  systems  effects  simulations  such  as  the  SLAD 
Scene  and  Countermeasure  Integration  for  Munition  Interaction  with 
Targets  (SCIMITAR)  program.  The  model  starts  with  a  known  given 
IR  scene  and  a  set  of  standard  meteorological  conditions  and  projects 
over  a  full  24-h  diurnal  cycle  or  longer.  It  is  required  that  the  standard 
meteorological  conditions  be  known  (from  either  a.posterioria  data  or 
from  other  models)  over  the  full  time  period  and  that  the  extended 
scene  area  be  reasonably  homogeneous.  The  model  addresses  the 
problem  in  two  parts:  (1)  to  predict  the  surface  temperature  (and 
hence  signature)  at  a  single  point  based  on  the  known  meteorology 
and  (2)  to  extrapolate  the  single  point  results  over  the  extended  scene 
area.  The  point  model  is  based  upon  the  equilibrium  "energy 
balance"  developed  by  Rachele  and  Tunick^  as  part  of  the  ARL 
Radiation  and  Energy  Balance  Field  Study  project.  The  model  also 
borrows,  in  part,  from  the  U.S.  Army  Corps  of  Engineers  Smart 
Weapons  Operability  Environment  (SWOE)  model^  for  characterizing 
the  extended  scene  area  in  terms  of  surface  thermal  and  radiative 
properties.  The  model  treats  effects  of  incoming  solar  and  sky 
radiation,  turbulent  heat  transport,  thermal  conduction  into  the  (soil) 
substrate,  and  evaporative  cooling.  Results  for  the  single  site 
experiment  agree  favorably  with  the  SLAD  Team  measurements  and 
extended  area  IR  scene  simulations  appear  reasonable.  Methods  for 
verifying  against  other  models  and  field  measurements  are  planned. 


1  Rachele,  H.  and  A.  Tunick,  "Energy  Balance  for  Imagery  and  Electro-magnetic 
Propagation,"  Journal  of  Applied  Meteorology,  33:8  (1994):  964-976. 

2  Welch,  J.P.,  "Smart  Weapons  Operability  Enhancement  Joint  Test  and  Evaluation 
Program:  Final  Report,"  U.S  Army  Corps  of  Engineers  SWOE  Report  94-1, 1994. 


IX 


1.  Introduction 


Most  of  iJie  existing  surface  "energy  balance"  models  used  in 
meteorological  applications  are  based,  for  the  most  part,  on  the 
following  rigorous  relationship  relating  temperature,  T,  to  the  energy 
flux  density,  F,  as 

PC,^^T^  =  -Vf(f,t)  (1) 

dt 

where  the  product,  pCp,  is  the  volumetric  heat  capacity  of  the 
medium,  t  is  time,  r  denotes  location  is  the  usual  Cartesian  coordinate 
designation,  emd  F  is  the  vector  heat  flux  density.  It  is  sometimes 
convenient  for  finite  element  modeling  to  express  eq  (1)  in  the 
following  equivalent  integral  form  as 

jpCp  dV  =  J(F  •  n)dA  (2) 

which  follows  directly  from  the  madiematical  Gauss  theorem  and 
applies  to  any  finite  volume  element  V  enclosed  by  the  corresponding 
surface  area  A  and  where  n  is  the  inward  directed  surface  area  normal 
unit  vector.  In  either  case  the  flux  density  vector  is  written  in  most 
general  form  as 

F(r,t)  =  R(r,t)  +  H(r,t)  +  S(r,t)  (3) 

where  the  three  terms  on  the  right  account,  respectively,  for  the  three 
major  heat  energy  transfer  mechanisms  of  radiation,  R,  convection,  H, 
and  conduction,  S. 

For  the  one-dimensional  case,  used  almost  universally  in  modeling 
the  surface  energy  balance  for  atmospheric  applications,  the  energy 
flux  vectors  are  all  directed  in  the  vertical,  in  which  case  upon 
inserting  eq  (3)  into  either  eq  (1)  or  eq  (2),  we  have  the  more 
convenient  scalar  form 

^  dT(z,t)  ,5R(z,t)  5H(z,t)  5S(z,t)^  ... 

which  is  valid  for  any  point  in  the  atmosphere  or  subsurface  as  long 
as  we  use  the  appropriate  values  for  the  specific  heat  capacity. 
Equation  (4)  as  written  does  not  include  any  direct  effects  of 
horizontal  advection  and  any  latent  heat  or  chemical  energy  storage 
terms,  but  is  otherwise  valid  for  any  point,  provided  we  used  the 
appropriate  values  for  the  thermal  properties  of  the  medium. 
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In  any  practical  model  it  is  convenient  to  separate  the  atmosphere- 
subsurface  system  into  the  three  regions:  the  air-atmosphere, 
subsurface  soil,  and  interface  layers  (fig.  1).  Applying  the  usual 
assumptions,  eq  (4)  is  expressed  in  each  of  the  three  subregions  as 


Atmosplieric  layer: 


(P  ^  p  )  ail 


dT(z.  t) 
dt 


aR(z,t)  ^  aH(z,t) 
5z  9z 


(5) 


Interface  region: 

,  dT(z,  t)  _aR(z,t)  ,  aH(z,t)  ,  5S(z,t) 

tP'-"plsfc  j.  ~  a  a  ai 

dt  OZ  OZ  OZ 


Subsurface  region: 


(pCp)sub 


dT(z,t) 

dt 


as(z,  t) 
dz 


(6) 

(7) 


where  we  have  explicitly  assumed  the  conduction  term  to  be 
negligible  in  the  air-atmosphere  layer  and  both  the  convection  and 
radiative  terms  to  be  negligible  in  the  soil  subsurface  layer  (as 
indicated  in  fig.  1).  The  interface  is  a  special  case  involving  all  three 
fluxes  that  are  dissipated  in  a  near  infinitestimal  layer. 


INTERFACE 


Figure  1.  Sketch  of  the  boundary  layer  model. 


The  expressions  for  the  atmospheric  and  subsurface  regions  (eq  (5) 
and  (7))  are  well  behaved  and  have  been  solved  in  various  degrees  of 
approximation  elsewhere  and  will  not  be  further  discussed  here 
except  when  needed  for  clarification.  Our  main  area  of  concentration 
is  on  the  interface  region  (eq  (6))  from  which  we  derive  the  surface 
temperature  that  ultimately  serves  as  a  botmdary  condition  for  the 
other  regions.  In  the  remainder  of  the  paper,  we  first  apply  the  above 
formalism  to  the  earth-air  interface  for  the  special  case  of  thermal 
equilibrium  using  the  general  approach  described  by  Rachele  and 
Tunick  [1]  and  then  examine  the  implications  of  the  approach  in  light 
of  the  more  rigorous  non-equiUbrium  formulation  proposed  by 
Sutherland.  [2]  We  then  extend  the  concept  to  include  an  entire  IR 
thermal  scene,  which  we  then  model  over  a  full  diurnal  cycle. 
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2.  Earth-Air  Surface  Formulation 


Concentrating  on  the  interface  region,  we  proceed  by  integrating  both 
sides  of  eq  (6)  over  an  arbitrarily  small  layer,  Az,  to  ultimately  arrive 
at  the  following  expression  for  the  interface  region  [2]: 

pC.^  =  R(0,t)  +  H(0,t)  +  S(Az,t) 

dt  (o) 

. «R(0,t)  +  H(0,t)  +  S(0,t) 

where  pCs  is  heat  capacity  per  unit  area  of  the  subsoil  surface.  As 
indicated  in  the  notation,  the  first  two  terms  on  the  right  represent  the 
radiative  and  convective  fluxes  as  evaluated  at  the  top  of  the  layer 
(z=0)  and  die  last  term  represents  the  conductive  flux  as  evaluated  at 
the  bottom  of  the  layer  (z=Az)  (see  fig.  2).  Here  and  throughout,  for 
the  interface  region  only,  we  treat  fluxes  directed  into  the  layer  as 
positive  and  fluxes  directed  out  of  the  layer  as  negative,  as  implied  in 
figure  2.  This  is  somewhat  different  ihan  the  usual  meteorological 
convention;  however,  it  is  convenient  and  even  necessary  in  order  to 
maintain  strict  consistency  with  physical  principles. 


Figure  2.  Expanded  sketch  of  the  earth-air  interface. 

The  next  step  is  to  determine  some  workable  mathematical 
formulations  for  the  various  fluxes  in  eq  (8)  with  a  desire  to  be  a 
simple  as  possible.  At  this  point  we  break  from  any  strictly  rigorous 
formulation  and  simply  accept  some  of  the  various  semi-empirical 
parameterization  schemes  developed  over  the  years  by  the  modeling 
community;  most  of  which  can  be  fotmd  in  contemporary  texts,  for 
example  those  by  Stull  [3]  and  Jacobs  [4]. 

We  begin  with  the  radiation  term,  R(z,t),  which  is  the  major  driver  of 
die  energy  balance  and  includes  both  a  downwelUng  shortwave 
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contribution  (~0.3-3.0um)  of  solar/ sky  origin  and  a  longwave 
contribution  (~3.0-200.0  um)  of  thermal  origin  which  includes  both 
the  downwelling  emissive  sky  contribution  and  an  upwelling 
emissive  surface  contribution,  all  evaluated  at  the  surface  and 
parameterized  as  follows: 

R(0,  t)  =  (0,  t)  +  R,,^  (0,  t)  +  R,,,  (0,  t) 

(9) 

. =  (nS.e-'''  +  nS,(l-e-'"))  +  e..,oTi-8.,.oT4 

where  the  shortwave  contribution  represented  by  the  first  two  terms 
(in  braces),  denoted  as  Rsun/  accounts  for  both  the  direct  solar  beam 
and  diffuse  scattering  over  the  full  upper  (sky)  hemisphere.  In  both 
terms  p  is  the  solar  zenith  angle  cosine,  is  the  solar  constant,  Sj  is  an 
empirical  parameter  for  diffuse  scattering,  and  r  is  the  atmospheric 
vertical  optical  depth  assumed  known  either  from  onsite 
measurements  or  from  other  existing  models.  Other  quantities  in 
eq  (9)  are  respectively  the  (known)  surface  emissivity,  Esfc ,  and  the  sky 
emissivity,  Csky  /  which  is  modeled  as  a  function  of  the  ambient  air 
temperature  and  relative  humidity  using  expressions  discussed  in 
appendix  A.  Note  in  eq  (9)  that  we  explicitly  define  the  upwelling 
surface  radiation,  Rsfc=-esfc(oTsfc)^,  as  negative  since  it  is  always 
directed  upward  and  away  from  the  surface.  Also  we  will  later  have 
occasion  to  approximate  the  "sky"  temperature  as  being  equal  to  the 
ambient  air  temperature  ( i.e.,  Tsky~Tair). 

The  convective  flux,  H(z,t),  also  called  the  sensible  heat  flux,  or 
"eddy"  heat  flux,  is  modeled  as  a  function  of  both  the  wind  speed, 
w(t),  and  the  air  surface  temperature  gradient  using  the  following 
semi-empirical  relationship: 

T  ftl-T  ftl 

H(0,  t)  =  K(w)  ^  (10) 

where  Tair(t)  is  the  (known)  air  temperature  and  K(w)  is  the  eddy 
diffusivity  which  is  a  function  of  the  measured  wind  speed,  w(f),  and 
can  be  inferred  from  boundary  layer  theory  or  other  empirical  models 
(cf.  eq  (18)  and  appendix  A).  Both  the  wind  speed  and  air  temperature 
are  assumed  to  be  measured  at  the  reference  height,  Azref,  which  for 
convenience,  we  take  to  be  one  meter  in  the  calculations  to  follow. 
Note  from  eq  (10)  that  when  Tair  is  greater  than  the  surface 
temperature,  Tsfc,  that  the  direction  of  the  flux  is  positive  (downward, 
toward  the  surface,  heating)  and  when  Tsfc  is  greater  than  Tair,  the 
direction  of  the  flux  is  negative  (upward,  away  from  the  surface, 
cooling),  all  consistent  with  the  stated  sign  convention. 
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The  subsurface,  or  soil,  conductive  heat  flux  density,  S(z,t),  is 
modeled  in  a  similar  way  using  the  following  simple  but  physically 
reasonable  representation: 

T  _  T  f  tl 

SiO,t)  =  c.X-^ - ^  (11) 

AZno, 

where  A  is  the  soil  thermal  conductivity  and  Tmax  is  the  "deep  soil" 
maximum  temperature,  which  is  assumed  a  known  constant;  and  cz  is 
an  adjustable  constant  determined  in  a  manner  described  later. 
Numerical  values  for  X  and  other  physical  parameters  as  used  here 
are  described  in  appendix  A.  Note  in  eq  (11)  that  when  Tsfc  is  less  than 
Tmax,  the  direction  of  the  flux  is  positive  (upward,  toward  the  surface, 
heating);  conversely  when  Tsfc  is  greater  than  Tmax,  the  direction  of  the 
flux  is  negative  (downward,  away  from  the  surface,  cooling),  all  again 
consistent  with  the  stated  sign  convention. 

Substituting  the  radiative  expression  into  eq  (8)  and  taking  care  to 
account  also  for  both  the  shortwave  reflection  and  longwave 
absorption  at  the  surface,  we  have 

pC,  ^  =  (1  -  r.,.)R_  +  +  H(0.t)  +  S(0,t)  (12) 

dt 

where  Ts/c  is  the  shortwave  surface  reflectivity  and  Ss/c  is  the  longwave 
surface  emissivity,  which  in  the  second  term  we  have  assumed  to  be 
equal  to  the  surface  absorptivity  in  accordance  with  the  well-known 
Kirchoff  approximation. 
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3.  Equilibrium  Surface  Fluxes 


There  is  some  controversy  regarding  fully  rigorous  time  dependent 
solutions  to  eq  (12)  because  the  details  in  the  infinitestimal  region  of 
the  interface  are  not  necessarily  well  understood.  However,  for  the 
equilibrium  case  we  have  dT /  dt=0  and  eq  (12)  can  be  solved  in  a 
number  of  ways  to  produce  what  we  will  call  here  the  equilibrium 
solution.  That  is,  setting  dT/dt=0  in  eq  (8)  we  have  immediately 

R,<,(t,T,,)  +  H,,(t,T,,)  +  S,,(t,T,,)  =  0  (13) 

where  we  have  inserted  the  subscript  "eq"  to  make  clear  that  the 
fluxes  represent  the  equihbrium  solutions.  We  now  seek  to  "balance" 
eq  (13)  by  finding  die  particular  value(s)  of  the  equilibrium 
temperature,  Teq(t),  that  produce  a  0  sum.  The  approach  taken  by 
Rachele  and  Tunick  [1]  is  based  on  an  iterative  scheme  to  find  the 
value  of  the  sensible  heat  flux  density,  Heq(t),  that  minimizes  the 
balance  in  eq  (13).  We  can,  in  principle,  use  these  results  to  express 
the  individual  equilibrium  fluxes  directly  in  terms  of  the  so- 
determined  equilibrium  surface  temperature: 

R.,  (t.T„)  =  (1  -  r^)R™  (t)  +  e.tR.1,  (t)  -  (t)  (14) 

T  -T  (t) 

Seq(t,Teq)  =  C,X  -  (15) 

^^not 

H.,(t.T^)  =  K(w)^*'’~^'^--  (16) 

where  all  parameters  remain  as  defined  in  earlier  expressions.  Some 
examples  of  how  the  method  rrught  work  using  measured  data  and 
the  above  formulation  are  shown  in  figure  3. 

The  plots  in  figure  3  are  divided  into  three  groups— upper,  middle, 
and  lower— and  are  based  on  data  gathered  from  a  set  of  special 
experiments  performed  at  the  White  Sands  Missile  Range  (WSMR)  by 
Anderson  and  Chenault  [5]  during  a  3-day  period  of  fair  weather 
conditions  in  the  late  fall  (Juhan  days  306,  307,  308).  As  shown  in  the 
upper  plots,  the  measured  data  included  air  temperature,  wind  speed, 
relative  humidity,  and  in  this  case,  the  measured  surface  temperature 
obtained  from  an  IR  (non-contact)  radiometer.  As  can  be  inferred 
from  inspection,  the  weather  during  the  3  days  was  reasonably  free  of 
any  adverse  events  and  showed  a  general  tendency  toward  increasing 
surface  and  air  temperatures.  The  only  direct  evidence  of  any  cloud 
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cover  is  the  "glitch"  in  surface  temperature  during  the  middle  of  the 
second  day  which  correlates  with  a  similar  event  in  the  solar 
measurement  (middle  plots).  The  wind  showed  a  slow  decline  during 
the  first  24  h  and  then  remained  near  a  steady  value  (about  2  m/  s)  for 
the  remainder  of  the  experiment.  The  relative  humidity  (shown  1/10 
scale  in  fig.  3)  consistently  increased  steadily  during  the  night  to  a 
maximum  of  around  60  percent  and  decreased  steadily  during  the 
day  to  a  minimum  of  around  20  percent,  roughly  synchronized  12  h 
out  of  phase  with  the  air  temperature. 


time(hrs) 

Figure  3.  DAY-24  meteorological  input  and  modeled  fluxes. 


For  these  experiments  the  shortwave  solar  input  was  measured 
directly  and  is  shown  in  the  middle  plots  of  figure  3  where  the 
previously  mentioned  "glitch"  (due  presumably  to  a  passing  cloud 
during  the  middle  of  the  second  day)  is  also  evident.  A  comparison  of 
the  measured  solar  data  with  the  model  of  eq  (8)  yielded  an  estimate 
of  t=0.215  for  the  atmospheric  optical  depth,  which  is  reasonable  for 
the  particular  site  and  season.  Although  it  is  hard  to  discern  from  the 
plots  the  trend  was  toward  a  steady,  but  small,  increase  in  solar  input 
with  time.  Shown  also  in  the  middle  plots  are  modeled  downwelling 
and  upwelling  longwave  radiative  fluxes  from  the  sky  and  surface  as 
derived  from  the  above  formulation  based  on  the  measured  air  and 
surface  temperatures.  Note  that  the  modeled  downwelling  longwave 
sky  radiation  is  nearly  constant  but  does  show  some  variation  due  to 


the  changing  relative  humidity.  The  upwelling  (surface)  radiation  is 
strongly  correlated  with  the  surface  temperature  and  is  generally 
larger  than  the  downwelling  (sky)  covmterpart  which  is  the  main 
mechanism  causing  the  surface  to  cool  at  night.  The  bold  trace  in  the 
middle  set  of  plots  is  the  net  radiative  flux  as  determined  from  the 
sum  of  all  three  radiative  contributions,  both  shortwave  and 
longwave.  It  is  this  net  radiative  flux  that  is  sometimes  referred  to  as 
the  "forcing"  term  that  drives  the  temperature  changes.  This  forcing  is 
generally  opposed  by  the  reaction  fluxes  from  the  air  and  subsurface 
(i.e.,  S  and  H)  which  we  have  shown  in  the  lower  set  of  plots  along 
with  the  net  contributions  formed  by  the  sum  of  the  two.  We  note  for 
this  particular  case  that  the  soil  heat  flux  was  generally  larger  (in 
magnitude)  than  the  air  heat  flux  at  night  but  less  during  the  day; 
however,  this  aspect  of  the  modeling  process  is  dependent  upon  the 
particular  parameterization  scheme  used  and  is  probably  the  greatest 
uncertainty  in  modeling  the  equilibrium  fluxes. 

To  elaborate  on  the  final  point  made  in  the  preceding  paragraph,  it  is 
clear  that  although  we  have  presumed  the  surface  and  subsurface  as  a 
simple  homogeneous  layer,  the  real  world  is  far  more  complicated 
and  actually  beyond  physical  description  in  a  rigorous  sense.  Our 
solution  to  dealing  with  this  important  practical  problem  lies  in 
determining  the  adjustable  parameter,  ca,  introduced  in  eq  (11).  In  our 
adaptation  we  use  the  a-posterioria  initial  measured  air-to-surface 
temperature  difference  to  adjust  die  value  of  ca  such  as  to  force  the 
equilibrium  to  zero.  That  is,  from  eq  (13)  and  (15)  we  have 


Ca  = 


Req  (to  >  Tn,ea  )  +  (t„ ,  T„^^  ) 

-FT  -T  (t  )1 

L'^max  *nieaV  o^J 


Az„ 


(17) 


where  Tmea  [=Tmea  (to)],  is  the  measured  initial  surface  temperature 
which  we  assume  to  be  at  least  approximately  equal  to  the 
equilibrium  temperature.  Once  ca  is  determined  from  the  initial 
conditions,  the  value  thus  obtained  is  maintained  constant 
throughout  the  full  3-day  period  for  which  the  model  was  run. 

The  particular  choice  made  here  in  forcing  the  equilibrium  balance  by 
adjusting  the  soil  heat  flux  is  not  the  only  option  available.  In  fact  the 
particular  value  obtained  in  this  manner  is  strongly  dependent  upon 
the  values  chosen  for  eddy  diffusivity,  K(w),  introduced  in  eq  (10) 
and  which  we  model  as  (see  also  appendix  A): 
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(18) 


K(w)  =  -K„( 


1  „  ,W(t)  +  W  1/2 


w, 


where  Ko  (=10.50  w-m-^-K-i)  is  the  zero  point  eddy  diffusivity  derived 
from  wind  tunnel  experiments  and  Ci  is  another  semi-empirical 
constant  to  account  for  the  effect  of  the  local  topography.  We  assume 
here  that  our  ci  parameter  can  be  associated  with  a  surface  roughness 
parameter  (ci~  Ln  (zruf/zo),  commonly  used  to  account  for  the  effect  of 
local  terrain  features  on  the  vertical  wind  profile.  In  the  initial  testing 
of  the  model,  we  found  a  value  of  Zmf  =20  (cm)  to  give  the  best  overall 
comparisons  with  the  a  posterioria  data  and  this  value  was  used 
throughout. 

One  way  of  actually  implementing  eq  (17)  to  determine  ci  if  the 
measured  initial  surface  temperature  is  not  known  is  shown  in 
figure  4. 


Figure  4.  DAY-24  effect  of  soil  thermal  properties. 
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The  middle  set  of  plots  in  figure  4  refer  to  the  soil  heat  flux,  S,  and  the 
upper  set  refers  to  the  sensible  heat  flux,  H  as  determined  from  eq  (15) 
and  (16)  assuming  the  actual  (measured)  surface  temperatures  to  be 
nearly  equal  to  the  equilibrium  temperatures.  The  various  subplots  in 
each  case  refer  to  guesses  of  the  initial  air-to-surface  temperature 
difference  [i.e.,  AT=Tair-Tsfc]  ranging  in  value  from  -10  to  0  °C  in  steps 
of  2.0  °C.  Since  the  initialization  begins  at  nighttime,  we  restrict  our 
guesses  to  the  inversion  case.  For  each  guess  we  calculate  a  new  value 
of  C2  based  on  eq  (17).  The  dark  solid  line  refers  to  the  "correct"  value 
based  on  the  measured  surface  temperature.  In  this  particular 
example  the  extreme  case  for  AT=-10  °C  for  the  soil  heat  flux  results  in 
a  nearly  flat  curve  near  zero  and  progresses  in  a  more  or  less 
monotonic  trend  upward  to  the  other  extreme  near  the  0.50  maximum 
for  AT=0.  For  this  particular  case  an  inversion  of  AT=-8  °C  represents 
tile  optimum  fit  which,  from  eq  (17),  gives  a  value  of  1.35  for  C2  which, 
in  turn,  corresponds  to  an  "effective"  thermal  conductivity  of  0.0041 
cal/mK-i.  However,  this  fitting  procedure  cannot  be  accomplished 
blindly  because  one  needs  also  to  look  at  the  effect  on  the  sensible 
heat  flux  which  in  this  particular  case  also  yielded  a  reasonable  curve, 
mainly  because  of  our  choice  for  Ci  in  eq  (18).  The  lowermost  plots  in 
figure  4  represent  the  corresponding  modeled  surface  temperatures 
and  are  discussed  later  in  section  4. 
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4.  Equilibrium  Surface  Temperature 


Although  it  is  simple  enough  to  calculate  the  surface  radiation  term 
directly,  it  is  convenient  for  computational  purposes  to  "linearize"  the 
energy  balance  expression  by  expanding  the  surface  radiation  term  in 
a  Taylor  series  and  approximating  as  follows  [4]; 

^sfc  ~  ^sfc^'l'sfc 

. =  +-^(T-TJ....  + higher  order  terms]  (19) 

oT 

. «  Esfc^lTlx  +  )] 

where,  following  the  usual  approach,  we  have  chosen  to  make  the 
expansion  about  Tmax,  the  (constant)  deep  soil  temperature.  Thus,  the 
"linearized"  version  of  eq  (12),  utilizing  eq  (19)  for  tiie  surface 
radiation  becomes 


pc,  ^  =  (1  -  )  R„  +  £.„R.„  -  e.,.o[Tl, + 4(Tl  (T-  T„ )] 

at 

+  W(T-T.„)  +  c,^(T-T„„) 


(20) 


Az 


ref 


Az, 


where  we  have  also  substituted  eq  (10)  and  (11)  for  the  last  two  terms. 
The  equilibrium  solution  can  be  found  by  setting  dT/  dt=0  in  which 
case  eq  (20)  can  be  solved  immediately  to  yield 


T„(t)  = 


(1  -  r.„)  R.„  (t)  +  etf.Ra,  (t)  +  +  K(w)  T,../  A  z„,  +  T.,,  (t)  /  Az„, 


Az 


(21) 


sub 


where  all  quantities  on  the  right  are  known  either  from  onsite 
measurement  or  from  the  various  parameterization  schemes 
docrunented  in  appendix  A.  Thus  eq  (21)  yields  an  immediate 
determination  of  the  equilibrium  surface  temperature  without  the 
need  for  extensive  reiteration. 


As  an  example,  we  show  the  equilibrium  solutions  in  the  lowermost 
set  plots  in  figure  4  (previous  section).  It  is  interesting  to  note  in 
figure  4  that  the  effect  of  the  soil  flux  parameter,  C2,  as  represented  by 
the  inversion  values,  is  to  lower  the  nighttime  minimum  and  at  the 
same  time  to  raise  the  daytime  maximum  which  is  a  consequence  of 
the  relationship  between  the  convective  and  conductive  fluxes.  In 
figure  4  we  also  show  (in  bold)  the  remotely  measured  surface 
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temperature  which  seems  to  match  best  with  an  initial  inversion  of 
-8  °C  which  was  actually  near  the  observed  value  of  -7.9  °C.  Note 
further  that  the  equilibrium  values  also  seem  to  accurately  model  the 
steady  increase  in  surface  temperature  as  observed  over  the  three  day 
period.  The  results  here  are  surprising  in  that  results  work  as  well  as 
they  do.  It  is  in  fact  something  of  a  paradox  that  the  time-independent 
solutions  do  even  a  reasonable  job  in  modeling  the  time  dependent 
scenario.  Some  insight  as  to  why  this  is  true  (or  not  true)  is  given  in 
the  next  section. 


5.  Time  Dependence 


d{(T(t) 


Although  it  is  tempting  to  assume  the  equilibrium  temperature  to  be 
the  same  as  the  actual  surface  temperature;  this  is  not  necessarily  true. 
To  find  the  actual  time-dependent  temperature  we  would  ordinarily 
need  to  start  over  with  the  time  dependent  version  of  eq  (8)  and 
proceed  directly  with  a  numerical  solution.  However,  there  is  more 
insight  to  be  gained  by  examining  the  theoretical  relationship 
between  the  actual  solutions  (i.e.,  dT/dt#0)  and  the  equilibrium 
solutions  (i.e.,  dT/dt=0). 

We  start  by  first  subtracting  the  equilibrium  solutions  represented  by 
eq  (13),  from  the  time  dependent  solutions  represented  by  eq  (8)  to 
obtain 


=  {R(t)  -  R,,(t)}  +  {H(t)  -  +  {S(t)  -  S,,(x)}  (22) 

t=T 

which  is,  of  course,  valid  since  the  equilibrium  solution  sums  to  zero. 

In  eq  (22),  we  have  introduced  the  dummy  time  parameter,  x,  to  avoid 
confusion  with  the  independent  time  variable,  t.  In  this  context  we 
need  to  think  of  the  equilibrium  values  at  any  time,  x,  as  constants 
based  upon  the  equilibrium  fluxes,  but  not  intrinsically  dependent  on 
time,  per  se.  Note  however  that  it  is  understood  that  the  derivative  is 
to  be  evaluated  at  t=x.  We  next  perform  the  term-by-term  subtractions 
using  the  basic  forms  for  the  various  fluxes.  This  operation  results  in 
some  cancellations  of  terms  and  we  ultimately  arrive  at  the  following 
final  expression: 


Te<,(x)} 


It 


pc, 


-{4aT, 


max  * 

Az 


+  W}{T(t)-T„(T)}  (23) 


Az 


ref 


where,  for  a  given  wind  speed,  each  of  the  three  terms  inside  the  first 
set  of  brackets  are  constants;  thus  eq  (23)  can  be  integrated 
immediately  to  yield 

AT(t)  =  (AT)„  exp-  {-1-  +  -L  4-  — }  t  (24) 

"^sfc  ^sub  ^air 


where  ATo  is  the  difference  between  tiie  actual  time  dependent 
temperature  and  the  equilibrium  temperature  [i.e.,  (AT)o=T(t)-Teq(  x)]. 
Thus,  eq  (24)  obviously  implies  an  exponential  relaxation  process 
with  three  distinct  time  scales  defined  as  follows: 
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pc, 

1  X  I  A  z 

^  sub  P  C  s 

1  _  K(w)/Az 

^  sir  pc. 


(25) 


where  the  first  accounts  for  radiative  processes,  the  second  for 
conductive  processes,  and  the  third  for  convective  processes. 

It  is  interesting  to  note  that  our  calculations  based  upon  the  model 
and  the  initial  measured  data  yielded  nominal  values  of  around  1.03, 
2.10,  and  1.51  min  for  each  of  the  respective  processes.  The  results 
here  indicate  that  the  time-dependent  solutions  will  be  different  in 
magnitude  from  the  equilibrium  values  and  there  should  be  some 
time  lag  due  to  the  various  inertial  processes  involved  and  this  does 
seem  to  be  evident  upon  close  examination  of  the  figure  4  data. 
Although  the  actual  correction  of  the  equilibrium  values  to  account 
for  the  time  dependence  is  beyond  the  scope  here,  we  are  pursuing 
the  matter  in  follow-on  work. 
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6.  Application  To  IR  Scene  Generation 


We  next  apply  the  energy  balance  concept  to  the  problem  of  IR  scene 
generation.  The  starting  point  is  some  initial  thermal  array  of  surface 
temperatures,  or  "map,"  such  as  that  shown  in  the  upper  graphic  of 
figure  5,  which  represents  a  "thermal"  image  for  the  nighttime  case 
near  local  midnight.  The  gray  scale  in  this  figure  is  arranged  such  that 
the  highest  temperatures  are  brighter  (i.e.,  white)  and  the  lower 
temperatures  are  darker  with  a  full  range  of  about  10  °C. 


basennap3Q8.dat 


2)  «  60  80  100  120  140  .  160  180 


Figure  5.  IR  scene  baseline  temperature  map  (upper)  and  soil  map  (lower). 


Note  in  figure  5  the  (cooler)  dirt  road  and  the  relatively  hotter 
vegetation  areas  that  are  characteristic  of  a  nighttime  scenario.  The 
particular  scene  in  figure  5  is  based  on  a  191  by  191  pixel  array 
corresponding  to  a  spatial  resolution  of  approximately  1  m  as  viewed 
from  directly  above.  The  data  were  collected  from  a  real  scene  which 
we  use  here  only  as  an  example  to  demonstrate  the  method. 


19 


The  bottom  array  of  figure  5  is  a  "thermal  characterization  map" 
which  represents  the  underlying  subsurface  thermal  properties  that 
were  generated  using  eq  (17)  and  pixel-by-pixel  air-to-surface 
temperature  differences  estimated  from  the  thermal  array.  In  the 
bottom  array  the  brighter  regions  correspond  to  higher  conductivity 
and  the  darker  regions  correspond  to  lower  conductivity.  In 
generating  the  soil  map  we  have,  of  course,  assumed  that  the  general 
weather  conditions  are  approximately  uniform  over  the  region  so  that 
we  can  use  the  single  "key  station"  data  to  model  the  various  fluxes  of 
the  energy  balance.  To  complete  the  example,  we  also  generated 
estimates  of  surface  reflectivity  and  emissivity  from  other  imagery  of 
the  region.  The  results  for  the  emissivity  and  reflectivity 
characterization  are  shown  in  figure.  6. 
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Figure  6.  IR  scene  emissivity  map  (upper)  and  reflectivity  map  (lower). 


The  characterizations  of  figure  6  are  based,  for  the  most  part,  on 
visual  band  imagery  used  to  discriminate  mainly  between  vegetated 
and  non-vegetated  areas  and  are  necessarily  crude.  In  both  cases  we 
used  three  levels  to  characterize  the  region  corresponding  to 
emissivity  values  of  0.90  for  bare  soil,  0.99  for  fully  vegetated,  and 
0.95  for  mixed  and  corresponding  reflectivity  values  of  0.20,  0.0,  and 


0.10,  respectively.  That  is,  note  in  the  upper  image  of  Fig.  6  that  the 
(shortwave)  reflectivity  in  the  vegetation  areas  is  generally  lower  than 
that  in  bare  soil  areas  as  we  might  expect  from  the  existing  database 
of  measurements.  Note  also  from  the  lower  image  that  the  (longwave) 
emissivity  in  the  vegetation  regions  is  generally  higher  than  that  in 
the  bare  soil  areas,  which  is  also  consistent  with  other  studies.  The 
characterization  here  is  necessarily  crude  and  is  intended  to  serve 
only  as  an  example. 

The  idea  now  is  to  use  the  onsite  point  data  from  the  key  location  as  a 
general  indicator  of  the  ambient  meteorology,  which  we  then  assume 
to  be  nearly  the  same  at  all  nearby  pixel  locations.  The  assumption 
then  is  ttiat  any  temperature  differences  must  be  due  to  the 
underlying  surface  properties.  Results  for  a  full  hour-by-hour  diurnal 
cycle  are  shown  in  figure.  7. 
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Figure  7.  IR  scene  diurnal  cycle. 


The  arrays  in  figure  7  begin  at  nearly  midnight  and  continue  through 
the  daytime  and  evening  until  midnight  the  next  day,  covering 
roughly  the  first  24  h  of  the  scenario  discussed  in  connection  with 
figures  3  and  4.  As  one  may  expect  from  examination  of  the  earlier 
results,  the  scenes  for  nighttime,  both  in  the  early  morning  and 
evening  periods,  appear  to  be  about  the  same,  indicating  slow 
changes,  which  is  consistent  with  the  meteorological  conditions  noted 
earlier.  However,  the  daytime  cases  are  more  interesting  and,  in  fact, 
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show  a  contrast  reversal  between  the  dirt  road  and  background 
between  about  0900  and  1500  hrs  including  the  neutral  event  at  1000 
and  1400  hrs  where  the  contrast  reduces  to  near  zero,  making  the  dirt 
road  nearly  undetectable.  Although  there  are  many  quantitative 
details  that  need  to  be  worked  out,  we  believe  that  the  results  here 
give  a  reasonable  qualitative  picture  of  reality  for  a  wide  variety  of 
applications.  Further  studies  are  in  progress  to  compare  results  with 
more  detailed  models. 


7.  Summary 


The  sxirface  temperature  results  for  the  single  site  example  compare 
favorably  with  measurements  provided,  indicating  that  we  made 
responsible  choices  for  the  various  empirical  parameters.  This 
statement  applies  both  to  the  various  fluxes  which  compare  favorably 
with  other  such  studies  in  the  literature  and  the  surface  temperature 
results  which  agree  with  data  from  our  own  experiments  at  WSMR. 
Further  development  and  experimentation  is  needed  to  establish  the 
accuracy  of  our  empirical  methods  for  modeling  the  underlying  soil 
properties  imder  various  meteorological  conditions.  Although 
appearing  realistic,  results  for  the  IR  scene  generation  need  to  be 
viewed  with  some  caution  as  we  have  not  yet  included  atmospheric 
effects  due  to  turbulence  and  aerosol-induced  obscuration  which  we 
plan  for  the  future. 
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Appendix  A.  Parameter  Definitions 


This  appendix  supplies  details  of  the  particular  parameterization 
schemes  used  to  mode  the  various  energy  fluxes  and  physical 
parameters,  including  subsurface  (soil)  thermal  properties,  eddy 
convection  parameters,  solar  shortwave  radiation,  sky  and  surface 
longwave  emission,  and  surface  moisture  evaporation.  There  is  often 
a  mix  of  units  used  throughout  the  literature,  and  it  is  worthwhile  to 
keep  the  following  conversions  and  definitions  in  mind: 

1  Calorie  =  4. 19002  Joules 
1  Langley  =  1  Calorie/minute  =  698.3366  watt/m^ 

1  mb  (millibar)  =  lOmPa  (milliPascal) 

Sg  =  Solar  constant  =  1.99  Langley  =  1353  watt/m^ 

=  latent  heat  of  vaporazition  =  2.45  *  10^  Joule  •  kg~' 

Rg  =  universal  gas  constant  =  287. 04 Joule  •  K~'  ■  kg~‘ 

=  gas  constant  for  water  vapor  =  46 1.5 Joule  •  K~‘  •  kg~‘ 

=  StefanBoltzmanConstant  =  5.67  *  10~^Watt  ■  m~^  •  K~^ 

Most  of  the  above  definitions  are  derived  from  the  standard  texts, 
such  as  those  by  StuU  [3]  and  Jacobs.  [4] 

We  start  with  the  subsurface  parameters  which  first  appear  in  eq  (11) 
and  which  we  model  as  a  function  of  soil  moisture  using  the 
following  empirically  determined  expressions: 

thermal  conductuvity : 

X,(x)  =  0.001 +  0.004  X (cal- cm"' •  sec"' •  K”’ ) 

specific  heat : 

pC(x)  =  0.27  + 1 .0  x*'" . (cal-  cm"'  -  K"' )  (Al) 

diffusivity : 

a{X)  =  ^ . (sec-' -cm-') 

pC(x) 

where  j  is  the  fractional  soil  moisture  content  and  p  is  the  soil  bulk 
density.  In  practice  the  actual  formulation  is  overridden  by  the  use  of 
the  reference  parameter  Azmax  and  the  adjustable  parameter  C2  in 
eq  (11);  however  the  expressions  do  give  us  some  type  of  check  with 
physical  entities.  In  the  Rachele-Tunick  [1]  model,  the  soil  heat  flux  is 
modeled  directly  as 
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(A2) 


S(t)  =  (T.„  -  T„  )^Sm[.^(t(hr)  - 1„ )  +  S. 

where  To,  to,  and  So  are  empirical  parameters  derived  from  the 
a.posterioria  data. 

The  eddy  diffusion  coefficient  first  introduced  in  eq  (10)  to  determine 
the  sensible  heat  flux  density  is  modeled  as  a  function  of  wind  speed 
using  the  following  empirically  determined  expression: 


=  =10.50 watt- m-'  -K'’  • 

c,  w, 


(A3) 


where  the  parameter  Ko  is  empirically  derived  from  wind  tunnel  data, 
w  (m/  s)  is  the  measured  wind  speed,  Wo  (=0.25  m/  s)  represents  the 
measurement  stall  speed,  wi  is  a  wind  scaling  parameter  taken  here  to 
be  unity,  mainly  to  keep  the  units  consistent,  and  Cj  is  an  empirical 
parameter  related  to  the  local  terrain  r 

oughness  as  ci=Ln(zruf/zo)  and  for  the  site  here  has  a  numerical  value 
of  2.95  [=Ln  20].  In  the  Rachele-Tunick  version  [1],  the  eddy  diffusion 
coefficient  is  not  used  directly  but  is  implicit  in  their  air-to-surface 
temperature  difference  formulation  expressed  as 


where  k  is  the  dimensionless  von  karman  constant,  T*  and  L  are 
scaling  parameters,  and  F^iz/L)  is  a  scaling  function  and  also  a 
dependent  on  the  local  roughness,  Zruf,  the  form  of  which  is 
determined  from  similarity  theory. 

The  sky  emissivity,  first  used  in  eq  (9)  is  modeled  as  a  function  of  the 
ambient  pressure,  temperature,  and  humidity  as: 

Sky  emissivity  (dimension  less): 

where  P(hPa)  is  the  ambient  atmospheric  pressure  and  eoap  is  the 
ambient  water  vapor  pressure: 


Ambient  vapor  pressure: 

_  RH(%) 

^vap- 


(1.831*10-’  )T,,,exp(0.06T,,3 ) . (hPa) 


(A6) 
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where  RH  (%)  is  the  ambient  relative  humidity  (in  percent)  and  Tats  is 
the  (absolute)  ambient  air  temperature. 

In  many  models  the  effect  of  surface  water  evaporation  is 
parameterized  as  an  equivalent  energy  flux  density,  usually 
designated  as  LE  and  as  being  proportional  to  the  sensible  heat  flux. 
This  suggests  a  new  parameterization  defining  the  combined  effect  of 
sensible  and  latent  heat,  H'  (=H+LE),  modeled  as 

H’(w,S„)  =  -%i^H(w)  (AT) 

X...,S„+T 

where  w  (m/ s)  is  the  ambient  wind  speed  and  H(w)  is  the  sensible 
heat  flux  density  (no  evaporation),  y  is  the  psychrometric  constant 
(=0.0004  K-i),  and  Xwet  is  a  dimensionless  empirical  parameter  varying 
from  a  value  of  0.0  for  a  completely  dry  surface  to  a  value  of  about 
1.25  for  a  nearly  saturated  surface.  In  eq  (A7)  the  quantity  Sec  is  the 
slope  of  the  saturation  vapor  curve  and  given  by 

S„  =  0.622  . (K'')  (A8) 

^gas'^abs 

where  Tabs  is  again  the  (absolute)  ambient  air  temperature,  Rgas  is  the 
universal  gas  constant,  Levap  is  the  latent  heat  of  vaporization  (see 
above  definitions  hst),  and  Qsat  is  the  saturation  specific  humidity 
given  by 

q„,=  0.622^ . (hPa)  (A9) 

where  Csat  is  given  by  eq  (A6)  and  P  is  again  the  ambient  pressure. 

We  next  turn  to  the  modeling  of  the  shortwave  radiation  originating 
from  the  direct  sun  rays  and  and  diffuse  sky  scattering.  The  solar 
zenith  Cosine,  p,  is  dependent  solely  upon  the  sxm-earth  geometry 
and  is  modeled  in  the  usual  manner  as 

Solar  zenith  angle  [ji=Cos(d)]: 

T  ,  -12  (AlO) 

p  =  Sin(0„,JSin(5,^)  +  Cos(e„,  JCos(53^)Cos(7t-^2!^) 

where  9„iat  is  the  measured  station  latitude  referenced  positive  north, 
Tsoiar  is  the  local  solar  time  (hours),  and  5sun  is  the  solar  declination 
angle  can  be  calculated  from  the  Julian  Day,  Jday,  as 
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Solar  declination  angle: 


T  -173 

Ssun  =  S^„Cos[360(-^^^— - )]... degrees  south 

3oj 


(All) 


where  6max  (=23.45°)  is  the  yearly  minimum  which  occurs  on  Julian 
day  173.  The  local  solar  time,  Tsoiar,  is  referenced  such  that  the  value 
Tsoiar=12  (hours)  corresponds  to  the  maximum  of  eq  (AlO)  which 
defines  the  station  local  noon.  The  relationship  between  solar  time 
and  Urdversal  Coordinated  Time  (UTC)  often  used  as  the  standard  in 
modeling  is 

Universal  Coordinated  Time: 

^  n  N  ^  ,  9ion  (degrees)  (A12) 

Tsoiar  (hours)  =  UTC(hours)  + - — - 


where  cpion  is  the  station  longitude  (0  to  ±  180°)  measured  as  positive 
progressing  eastward,  and  negative  progressing  westward  from  the 
prime  meridian  at  Greenwich,  UK.  The  solar  constant.  So,  is  a  function 
of  the  earth-sun  distance  which  undergoes  some  variation  throughout 
the  year  and  is  modeled  as 

S<.(J  )  =  S.  { 1  +  0.338  Cos[^^%y2)])  ....(watt,  nf’ )  (A13) 

ioj 

where  So  is  the  mean  solar  constant  (1353  watt-  m’^ ). 
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